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Abstract.  The  Bareiss  (or  Schur)  and  Levinson  algorithms  are  the  mo^t  popular  algorithms  for 
solving  linear  systems  with  dense  n  x  n  Toeplitz  coefficient  matrix  in  O(n')  arithmetic  operations. 
Both  algorithms  have  been  generalised  to  solve  linear  systems  whose  n  X  n  coefficient  matrices  A  are 
not  necessarily  Toeplitz  (in  0(n3)  operations).  We  show  in  this  paper  that  the  generalised  Levin¬ 
son  algorithm  is  a  direct  consequence  of  the  generalised  Bareiss  algorithm,  thereby  considerably 
simplifying  its  presentation  in  comparison  to  previous  work. 
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1.  Introduction 

The  solution  of  linear  systems  of  equations  Ax  =  b  with  real  dense,  non-singular  coefficient 
matrix  A  is  usually  accomplished  by  first  determining  a  triangular  factorisation  of  A  [9].  There  are 
two  types  of  triangular  factorisations:  lower-upper  and  upper-lower.  A  lower-upper  factorisation 
is  given  by  A  =  LU ,  where  L  is  a  unit  lower  triangular  and  U  an  upper  triangular  matrix,  while 
an  upper-lower  factorisation  provides  A  =  UL,  where  U  is  unit  upper  triangular  and  L  lower 
triangular  (a  unit  triangular  matrix  is  a  triangular  matrix  with  ones  on  the  main  diagonal).  For 
non-singular  A  the  factorisations  are  unique.  The  number  of  arithmetic  operations  required  to 
factor  a  n  x  n  matrix  is  0(n3). 

In  Section  2  we  give  a  simple,  structural  view  of  these  triangular  factorisations,  which  is 
intended  to  lead  to  and  facilitate  the  explanation  of  the  generalised  Bareiss  algorithm  in  Section  3. 
The  generalised  Bareiss  algorithm  [5]  computes  a  matrix  Q  and  the  two  triangular  matrices  U  and 
L  such  that 
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As  a  product  of  2  x  2  elementary  transformations,  Q  is  of  the  form 


where  L\  +  Li  is  unit  lower  triangular  and  U\  -f  Ui  unit  upper  triangular.  The  uniqueness  of  the 
factorisations  implies  Z_1  =  L\  +  Li  and  U\  +  Ui  =  U~l.  Thus,  the  generalised  Bareiss  algorithm 
computes  explicitly  the  triangular  factors  U  and  L  of  A.  The  unit  triangular  factors  L~x  and  U~l  of 
A-1  are  obtained  implicitly:  in  factored  form  as  the  product  Q  of  2  x  2  elementary  transformations. 

In  Section  4  the  generalised  Bareiss  algorithm  is  scaled  to  yield  a  symmetric  factorisation  for 
symmetric  positive-definite  matrices  A,  the  so-called  Hyperbolic  Cholesky  algorithm  [4]. 

Upon  multiplication  of  the  components  in  equation  (*)  by  A-1  one  obtains 

«(!)-(£)• 

where  I  is  the  identity  matrix.  This  is  the  generalised  Levinson  algorithm  [6]  of  Section  5;  it 
explicitly  computes  the  matrix  Q  and  the  unit  triangular  factors  L-1  and  U~l  of  A-1.  The  two 
triangular  factors  of  A  are  available  implicitly,  they  may  be  determined  from  U  —  L~lA  and 
L  =  U~lA. 

We  also  discuss  how  the  matrix  Q  can  be  used  to  accomplish  forward  elimination  and  back- 
substitution  in  order  to  determine  the  solution  x  to  the  linear  system  Ax  =  6,  when  A  is  symmetric 
positive-definite. 

Section  6  concludes  the  paper  with  a  combination  of  the  Bareiss  and  Levinson  algorithms  for 
the  fast  parallel  solution  of  linear  systems  with  persymmetric  coefficient  matrices.  For  the  special 
case  of  Toeplitz  matrices,  it  seems  to  be  identical  to  an  algorithm  from  [8],  but  we  provide  a  much 
simpler  description. 

In  each  section,  we  discuss  how  the  algorithms  take  advantage  of  the  situation  where  the 
coefficient  matrix  A  is  a  Toeplitz  matrix.  A  n  x  n  matrix  T  is  a  Toeplitz  matrix  if  its  elements  are 
constant  along  the  diagonals: 
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Linear  systems  of  equations  Tx  =  6  with  Toeplitz  coefficient  matrix  T  arise,  for  instance,  from 
time  series  analysis,  Pade  approximations,  and  discretisations  of  partial  differential  and  integral 
equations  [3].  Because  a  nxn  Toeplitz  matrix  contains  at  most  2n—  1  different  elements,  Toeplitz 
linear  systems  can  be  solved  with  0(n 2)  or  fewer  arithmetic  operations.  In  this  case,  the  generalised 
Bareiss  and  Levinson  algorithms  reduce  to  the  classical  Schur  [2,  13]  and  Levinson  [12]  algorithms, 
which  perform  the  solution  ofnxn  Toeplitz  systems  in  0(n2)  operations. 

The  purpose  of  this  paper  is  to  provide  simple  algebraic  explanations  for  the  algorithms.  We 
hope  that  our  strategy  of  first  presenting  the  methods  for  non-Toeplitz  matrices  followed  by  the 
subsequent  restriction  to  Toeplitz  matrices  gives  a  much  clearer  idea  of  the  underlying  structure  of 
the  algorithms.  The  numerical  issues  associated  with  the  methods,  which  we  ignored  in  this  paper, 
are  discussed  in  [3,  9].  We  only  consider  here  the  case  of  real  matrices,  but  the  methods  generalise 
readily  to  the  complex  case.  Reference  [9]  contains  all  the  relevant  facts  about  matrix  theory. 

2.  Solution  of  Linear  Systems  based  on  Gaussian  Elimination 

Let  A  be  a  real  non-singular  coefficient  matrix  and  b  the  right-hand  side  vector  of  the  system 
of  simultaneous  linear  equations  Ax  =  b.  When  A  is  a  dense  matrix  such  systems  are  usually  solved 
by  computing  a  triangular  factorisation  of  A. 

2.1.  General  Matrices 

There  are  essentially  two  types  of  triangular  factorisations.  If  all  leading  principal  submatrices 
of  A  are  non-singular  then  the  factoriation  A  —  LU  exists  and  is  unique,  where  L  is  a  unit  lower 
triangular  matrix  (i.e.  a  triangular  matrix  with  ones  on  the  diagonal)  and  U  is  an  upper  triangular 
matrix  with  non-zero  diagonal  elements.  The  linear  system  Ax  =  b  can  then  be  solved  in  three 
steps 

1.  Factorisation  A  =  LU 

2.  Forward  Elimination  Lc  =  b 

3.  Backsubstitution  Ux  =  c. 

Alternatively,  if  all  trailing  principal  submatrices  of  A  are  non-singular,  then  the  factorisation 
A  -  UL  exists  and  is  unique,  where  U  is  a  unit  upper  triangular  matrix  and  L  is  a  lower  triangular 
matrix  with  non-zero  diagonal  elements,  so  Ax  =  b  may  be  solved  by  computing 

A  =  UL,  Ud  =  6,  Lx  =  d. 

Applying  Gaussian  Elimination  (without  pivoting)  to  factor  the  nxn  matrix  A  requires  ^n3+0(n2) 
arithmetic  operations,  and  performing  the  two  subsequent  steps  by  triangular  system  solution 
requires  another  0(n2)  operations. 

In  many  applications  the  matrix  A  satisfies  a  uger  condition  than  the  two  above:  it  is 
symmetric  positive-definite  (and  so  are  all  its  contiguous  principal  submatrices).  In  this  case,  U  = 
DiLt  and  L  =  DyUT  and  the  two  symmetric  factorisations  A  =  LDi,Lt  and  A  =  U D\jUT  exist, 
where  and  Z?u  are  diagonal  matrices  with  positive  elements  on  the  main  diagonal.  Symmetry 
saves  half  of  the  work  in  the  factorisation  step,  so  only  ^n3  +  0(n2)  arithmetic  operations  are 
necessary. 

2.2.  Toeplitz  Matrices 

The  Levinson  algorithm  exploits  the  persymmetry  of  Toeplitz  matrices  T:  Tt  =  JTJ ,  where 
J  is  the  permutation  matrix,  also  called  ‘exchange  matrix’,  with  ones  on  the  anti-diagonal  [9]. 
Obviously,  if  T  is  persymmetric,  so  is  its  inverse. 

Persymmetry  implies  that  the  triangular  factors  from  the  two  types  of  factorisations  are  per¬ 
mutations  of  each  other.  In  particular,  let  U  =  where  Z?l  is  a  diagonal  matrix  and  Ui,  a 
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Figure  1:  Evolving  Non-Zero  Structure  of  a  4  X  4  Matrix 
During  the  Computation  of  A  =  LU. 

unit  upper  triangular  matrix,  and  similarly  L  =  D\}L\j ,  where  Dyj  is  diagonal  and  L\j  unit  lower 
triangular.  Then 

T  =  LDlUl  =  UDvLu,  and  T  =  JTtJ  =  ( JL%J)(JDVJ)(JUTJ ), 

implying  that  JLj^J  =  L,  JD\jJ  =  Di,  and  JUTJ  =  Ui,. 

If  T  is  also  symmetric,  i.e.  t{  =  t_,  for  1  <  i  <  n  —  1,  then 

JTJ  =  T  =  LDlLt  =  UD\jUt, 


so  L  —  JU  J.  In  particular,  the  last  row  of  T  is  the  reverse  of  its  first  row,  and  its  last  column  the 
reverse  of  its  first  column. 

The  solution  methods  outlined  above  are  not  able  to  exploit  any  Toeplitz  structure.  That  is, 
even  though  the  n  x  n  matrix  T  is  described  by  at  most  2n  -  1  parameters,  the  operation  count  for 
the  linear  system  solution  is  still  0(n3).  In  order  to  explain  algorithms  that  solve  Toeplitz  systems 
in  0(n2)  arithmetic  operations,  it  is  helpful  to  first  consider  how  Gaussian  elimination  changes  the 
zero  structure  of  a  general  matrix  during  the  factorisation  process. 

2.3.  Structural  View  of  the  Factorisations 

Gaussian  elimination  (without  pivoting)  accomplishes  the  factorisation  A  =  LU  by  computing 
successive  columns  of  the  unit  lower  triangular  matrix  L  and  successive  rows  of  the  upper  triangular 
matrix  0  by  zeroing  out  successive  columns  in  A.  The  zero-structure  of  the  matrices  occuring  in 
the  evolution  from  the  original  matrix  A  to  the  final  upper  triangular  U  are  depicted  in  the  4x4 
example  of  Figure  1.  There,  i  represents  an  element  which  is  generally  non-zero,  u  an  element 
of  the  final  matrix  (/,  a  blank  a  zero  element  and  ®  an  element  doomed  for  elimination  in  the 
current  step.  The  partitioning  distinguishes  those  elements  in  the  upper  part  of  the  matrix  that 
have  ceased  to  participate  in  computation.  In  each  step,  an  appropriate  multiple  a  of  the  latest 
row  of  (/,  which  is  the  one  just  underneath  the  horizontal  partitioning  line,  is  subtracted  from  each 
row  beneath  it  in  order  to  remove  the  elements  in  the  next  column.  So,  all  operations  are  of  the 
form 

lower  row  =  lower  row  —  a  *  latest  row  of  U. 

In  this  process  no  previously  introduced  zeros  are  destroyed,  because  the  latest  row  of  tJ  has  the 
same  zero  structure  as  all  the  lower  rows  to  which  it  is  added.  Note  that  in  each  step  all  the  rows 
can  be  modified  simultaneously  and  independently. 

The  computation  of  the  factorisation  A  —  UL  proceeds  in  a  similar  manner. 

3.  The  Generalised  Bareiss  Algorithm 

In  [2]  Bareiss  proposed  an  algorithm  for  solving  square,  non-symmetric  Toeplitz  systems.  The 
process  of  factoring  the  matrix  in  Bareiss  algorithm  is  identical  to  the  Schur  algorithm,  which  is 
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Figure  2:  Evolving  Non-Zero  Structure  of  a  4  x  4  Matrix 
During  the  Generalised  Bareiss  Algorithm. 


based  on  [13];  in  addition,  Bareiss  realised  that  forward  elimination  could  also  be  accomplished  by 
recursions  similar  to  those  used  for  the  factorisation.  In  [5]  it  was  shown  how  to  generalise  this 
algorithm  to  non-Toeplitz  matrices.  In  contrast  to  Gaussian  elimination,  which  computes  either 
A  =  LU  or  A  =  UL ,  the  generalised  Bareiss  algorithm  computes  factors  of  both  factorisations 
A  =  LU  and  A  =  UL  by  working  with  two  copies  of  the  matrix  A: 


Q 


The  next  section  explains  how  to  reduce  two  copies  of  the  original  matrix  A  to  an  upper  triangular 
matrix  and  a  lower  triangular  matrix,  and  the  one  following  argues  that  these  triangular  matrices 
are  (in  exact  arithmetic)  identical  to  the  ones  obtained  from  Gaussian  elimination. 

3.1.  Structural  View  of  the  Factorisation 

In  contrast  to  Gaussian  elimination,  which  eliminates  columns,  the  generalised  Bareiss  algo¬ 
rithm  eliminates  diagonals,  which  turns  out  to  be  crucial  for  the  exploitation  of  Toeplitz  properties. 
This  is  shown  in  the  4x4  example  of  Figure  2;  there  an  element  of  the  final  matrix  L  is  denoted 
by  /,  and  the  letters  x,  y,  z  denote  elements  that  are  generally  non-zero.  The  upper  copy  of  A  is 
transformed  to  upper  triangular  form  and  the  lower  copy  to  lower  triangular  form.  This  is  accom¬ 
plished  by  ‘rotating’  in  step  k  rows  i  +  k  in  the  upper  matrix  with  with  rows  i  in  the  lower  matrix 
(in  Figure  2  those  rows  are  made  up  of  identical  letters  x ,  y  or  z)  so  as  to  eliminate  one  element  in 
each  of  them,  namely  (i  +  k,  i)  and  (t,  i  +  k ): 

(  row  i  +  k  in  upper  matrix^  _  /  1  -Oi.i+fc  "\  (  row  i  +  k  in  upper  matrix 

\  row  i  in  lower  matrix  )  ^  —(3t+k,i  1  /  \  row  i  in  lower  matrix 

The  two  rows  are  selected  to  have  the  same  zero  structure,  so  already  introduced  zeros  are  preserved. 
Note  that  in  each  step  all  pairs  of  rows  can  be  modified  simultaneously  and  independently.  The 
transformation  matrix  Q  is  thus  made  up  as  a  product  of  these  2x2  ‘rotations’. 

We  can  now  express  the  generalised  Bareiss  algorithm  formally.  As  a  matter  of  simplicity,  we 
just  refer  to  entire  rows  of  the  matrix  rather  than  individual  elements  and  do  not  show  how  one  can 
save  operations  by  taking  advantage  of  the  increasing  number  of  zero  entries  in  the  matrix.  Below, 
the  vectors  a-0*  are  initialised  to  be  the  t'th  row  of  the  matrix  A.  From  now  on,  quantities  carrying  a 
positive  superscript  will  be  associated  with  the  upper  matrix,  and  those  with  a  negative  superscript 
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will  be  associated  with  the  lower  matrix.  In  the  A:th  step,  the  kth  subdiagonal  of  the  upper  matrix 
and  the  &th  superdiagonal  of  the  lower  matrix  are  removed;  that  is,  elements  (t  +  k,  i )  in  the  upper 
matrix  and  elements  (i,  i  +  k)  in  the  lower  matrix  are  removed  by  appropriately  combining  the 
(i  +  fc)th  row  a-* fc1)  of  the  upper  matrix  with  the  ith  row  of  the  lower  matrix. 

Factorisation  in  the  Generalised  Bareiss  Algorithm 

1  <  *  <  n,  a-0^  =  a, 

1  <  k  <  n  —  1,  l  <  i  <  n  —  k, 


( 4 i+fc  \_(  i  \  ( 45 

Us-fcV"v-^.«  1  A  «!" 


(*-l) 
i+k 
-k+ 1) 


Now  the  fcth  rows  of  U  and  L  are  respectively  given  by  a[.fc  ^  and  n+*). 

3.2.  Correctness  of  the  Factorisation 

It  remains  to  explain  why  indeed  the  generalised  Bareiss  algorithm  computes  triangular  fac¬ 
torisations  of  the  matrix  A. 

As  illustrated  in  the  previous  section,  the  transformation  matrix  Q  that  reduces  the  two  copies 
of  a  A  to  triangular  matrices  consists  of  products  of  2  x  2  ‘rotations’,  and  ‘rotations’  belonging  to 
the  same  step  are  disjoint.  The  matrix  Q  for  the  4x4  example  in  Figure  2  is  of  the  form 


1 

1 

1 

1 

X 

X 

1 

l 

1 

1 

The  right-most  of  the  three  matrices  eliminates  the  first  pair  of  off-diagonals.  The  x’s  in  the  first 
subdiagonal  of  its  top  right  part  represent  the  multipliers  for  removing  the  first  subdiagonal  in  the 
upper  matrix;  similarly  the  x’s  in  the  first  superdiagonal  of  the  bottom  left  part  are  the  multipliers 
for  removing  the  first  superdiagonal  in  the  lower  matrix. 

Because  subdiagonals  are  removed  in  the  upper  matrix  and  superdiagonals  are  removed  in 
the  lower  matrix,  the  transformation  matrices  applied  at  each  step  consist  of  two  lower  triangular 
matrices  in  the  top  half  and  two  upper  triangular  matrices  in  the  bottom  half  of  the  matrix.  Due 
to  the  particular  zero  structure  of  these  triangular  matrices  the  product  Q  of  these  transformation 
matrices  is  shown  in  [5]  to  have  the  form 


_(LX  L2\ 
-\UX  Ut) 
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If  /I  is  a  n  x  n  matrix,  then  the  2  n  x  2n  matrix  Q  consists  of  four  nxn  triangular  matrices,  whereby 
L\  is  unit  lower  triangular,  L2  is  strictly  lower  triangular,  Ux  is  strictly  upper  triangular,  and  U2  is 
unit  upper  triangular.  But  Lx  +  Z,2  is  a  unit  lower  triangular  matrix,  and  so  is  its  inverse;  similarly 
U\  +  U2  and  its  inverse  are  unit  upper  triangular  matrices.  Moreover,  in 

o  ( +  L2)a\ 
v  w "  \(Ux  +  u2)a) 

(L 1  +  Ij2)A  is  upper  triangular,  and  {Ux  +  Ui)A  is  lower  triangular.  Due  to  the  uniqueness  of  the 
triangular  factorisations  we  must  then  have  L  =  (Lx  +  I2)-1,  U  =  (Ux  +  U2)~x,  (Lx  +  L2)A  =  U 
and  (Ux  +  U2)A  =  L .  Consequently, 

Thus,  the  generalised  Bareiss  determines  explicitly  the  triangular  factors  U  and  L  of  A  and 
implicitly  (as  the  product  of  the  ‘rotations’  making  up  Q )  the  unit  triangular  factors  L~x  and  U~x 
of  /l-1. 

In  [5]  it  is  proved  that  the  generalised  Bareiss  algorithm  does  not  break  down  if  all  contiguous 
principal  submatrices  of  A  are  non-singular.  This  is  a  more  stringent  condition  than  needed  for  the 
independent  computations  of  A  =  LU  and  A  =  UL  by  Gaussian  elimination,  which  only  require 
all  leading  or  all  trailing  principal  submatrices  of  A  to  be  non-singular. 

3.3.  Forward  Elimination  and  Backsubstitution 

Instead  of  solving  a  triangular  system  the  transformation  matrix  Q  may  be  used  to  perform 
forward  elimination,  as  Bareiss  [2]  already  realised  in  the  Toeplitz  case.  From  the  last  equation  in 
the  previous  section  we  see  that  the  application  of  Q  amounts  to  a  premultiplication  by  L~x  and 
U~x,  thus  Q  can  also  be  applied  to  the  right-hand  side  b  in  order  to  effect  forward  elimination 

Backsubstitution  can  now  be  performed  by  solving  either  Ux  =  c  or  Lx  =  d  in  \  +  C(n) 
arithmetic  operations.  Employing  a  trick  from  Bareiss  [2],  one  could  alternatively  determine  the 
upper  half  of  x  from  Ux  =  c  and  the  lower  half  from  Lx  =  d,  which  would  require  only  ^  +  0(n ) 
operations. 

If  A  is  symmetric  positive-definite,  the  application  of  QT  can  replace  the  triangular  system 
solution  for  backsubstitution.  To  this  end,  let  I  denote  the  identity  matrix  of  the  same  size  as  A, 
and 

y  =  (I  I)Qt(^D^  D-1^(<j)=(Lx  +  L2)TDZxc  +  {U1  +  U2)TDuld 

=  L-TD£lc+U-TDvXd 
because  c  =  L~xb  and  d  =  U~xb,  so 

y  =  L-TD^L~xb  +  U-TDvXU~xb  =  A~xb+  A~xb  =  2^_16  =  2x. 

In  order  to  obtain  x,  rather  than  2x,  observe  that  c  and  d  each  give  rise  to  one  x,  and  applying  the 
transformations  instead  to  ( Ac  (1  -  A)d)r,  where  A  is  any  real  number,  yields 

y  =  (I  /)Qr  1)  ((1-CA)d)  =  A^i  +  U-A)/!-1^  Ax  +  (1-A)x  =  x. 


6 


* 


<0  <1  tj  <3 

t-1  to  tl  f? 

t-2  t-1  to  tl 

t— 3  t— 2  t_i  to 

to  ti  tj  h 

l0  l2 

t •  t t  *• 

—2  f0  £1 

V  t*  t* 

1  —3  —2  l0 

to  tj  tj  t3 

tf  t '  t‘ 

l0  l2 

*0  ll 

t"  t " 

l-3  l0 

— 

to  ti  t2  (3 

r0  ll  l2 

fit  tn 
c0  M 

l0 

30  *1  32  33 

S—  l  So  3l  32 

3- 2  »-l  30  31 

.  3-3  3-2  S_1  30  . 

9 0  92  53 

S-1  50  s2 

a'  a '  a1 

3 -2  -1  90 

4-3  3-2  4-1  4q 

a"  a" 

30  *3 

a-2  *-l  *0 

5-3  3-2  9-l  a0 

a"f 

90 

3-1  *0 

*'-2  *‘-1  *0 

3- 3  3-2  J-l  90 

Figure  3:  Evolving  Structure  of  a  4  X  4  Toeplitz  Matrix 
During  Bareiss  Algorithm;  Initially  s,  =  t,. 


In  particular,  A  may  be  set  to  either  zero  or  one  so  as  to  require  only  one  of  the  right-hand  sides  c 
or  d  as  input. 

For  symmetric  positive-definite  matrices  A  the  generalised  Bareiss  algorithm  permits  the  so¬ 
lution  of  the  linear  system  Ax  =  6  with  the  following  three  steps: 


1. 


Determine  Q  such  that  Q 


2.  Determine 


Determine  i  =  (  /  I)QT (  V  )  ((1  ^A)i)  • 

Because  the  generalised  Bareiss  algorithm  determines  a  non-symmetric  factorisation  <*,,,+*  / 
Pi+k,i,  and  it  can  generally  not  take  advantage  of  the  symmetry  of  a  matrix  in  order  to  reduce  the 
number  of  arithmetic  operations  -  except  if  the  matrix  is  a  symmetric  Toeplitz  matrix:  then  the 
operation  count  decreases  by  fifty  percent  compared  to  the  non-symmetric  Toeplitz  case,  as  shown 
in  the  next  section. 


3.4.  Specialisation  to  Toeplitz  Matrices 

If  the  matrix  T  to  be  factored  is  a  Toeplitz  matrix,  then  the  generalised  Bareiss  algorithm  is 
identical  to  the  original  Bareiss  algorithm  [2],  which  exploits  the  structure  of  a  Toeplitz  matrix. 
Figure  3  displays  the  evolution  of  the  structure  of  the  4x4  Toeplitz  matrix 


/  to 

h 

1 2 

h\ 

t- 1 

to 

1 1 

t2 

t-2 

t-1 

to 

t\ 

'*-3 

t-2 

t-l 

to  / 

during  the  course  of  Bareiss  algorithm.  Although  the  triangular  factors  of  a  Toeplitz  matrix  are 
generally  not  Toeplitz,  the  ‘working  part’  of  the  matrix  in  Bareiss  algorithm  is  Toeplitz  as  shown 
in  Figure  3  and  below,  because  identical  elements  along  each  diagonal  result  in  identical  Qi+k,i  and 
for  each  step;  this  makes  it  possible  to  exploit  the  Toeplitz  structure. 

In  particular,  denote  the  elements  of  T  by  =  tj_,,  1  <  i,j  <  n.  Rows  a\  and  an  constitute 
the  respective  first  row  of  U  and  last  row  of  L.  Since  T  is  Toeplitz,  at  the  beginning  of  the  first  step 
all  the  elements  necessary  for  further  computation  are  contained  in  rows  a?  and  an  of  the  upper 
matrix  and  rows  ai  and  an_i  of  the  lower  matrix.  The  multipliers  are 


_  h 

ti,i  Iq 


=  <*1, 


Pi+l,i 


ti.i+l  _  (-l  _ 
*,+!., +1  t0 


1  <  *  <  n  —  1. 
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These  multipliers  are  used  in  the  rotations,  as  above,  to  construct  rows  aJ+j  of  the  upper  matrix 
and  rows  a[-1^  of  the  lower  matrix,  1  <  t  <  n  -  1.  Hence  rows  3,...,n  of  the  upper  matrix  and 
rows  1,. .  .,ti  —  2  of  the  lower  matrix  retain  their  Toeplitz  structure,  and  rows  a^1*,  al1*,  and 

encompass  all  necessary  information  for  further  computation,  see  Figure  3. 

In  general  at  the  beginning  of  step  fc,  rows  k  +  1,.  ..,n  of  the  upper  matrix  as  well  as  rows 
l,...,n  —  k  of  the  lower  matrix  have  Toeplitz  structure.  If  we  represent  the  upper  and  lower 
triangular  parts  of  the  two  matrices  by  separate  vectors 


^  =  (o 

6(rfc+1)  =  (  o 


ak+l,k+l 

J-k+D 

al,k+l 


“jfc+t.n/’  °n  V“n,  1 

al,n  /  ’  °n-k  ~  V  an-k,l 


,(*-1) 

*n,n—  k 


0 


0) 


(-*4-1) 

n— A:,n— it 


o), 


then  and  ^  contain  all  the  distinct  elements  in  the  respective  upper  and  lower  triangular 

parts  of  the  upper  matrix,  while  and  hj-*r+1^  contain  all  distinct  elements  in  the  respec¬ 

tive  lower  and  upper  triangular  parts  of  the  lower  matrix,  see  Figure  3.  Thus,  the  following  two 
‘rotations’  embody  all  the  distinct  computations  on  the  Toeplitz  matrix  during  step  k: 


The  transition  to  the  next  step  is  accomplished  by  considering  the  next  lower  row  in  the  upper 
matrix  and  the  next  higher  row  in  the  lower  matrix: 


C**i2i*  bn-k+i=bn-kZT,  where  Z  = 


/0  1 


V 


0  1 


0/ 


This  leads  to  the  following  form  of  the  factorisation  process  in  Bareiss  Algorithm,  which  is  also 
called  ‘Schur  algorithm’  [13].  Below,  0*  denotes  a  row  vector  consisting  of  k  zeros,  and  the  rows 
are  represented  by  their  individual  elements. 

Schur  Algorithm:  Factorisation  Part  of  Bareiss  Algorithm 
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Now  the  fcth  rows  of  U  and  L  are  respectively  given  by 


(ofc  t 


(*)  tW 


t{k)  ) 

•  •  •  ln-k- 1  > 


and 


\ t-n+fc+l 


4"fc)  0^ )  • 


The  recursions  for  forward  elimination  and  backsubstitution  are  similar  to  those  in  the  factorisation; 
see  also  [4]. 

If  the  Tf'  ^itz  matrix  T  is  symmetric  positive-definite,  then  aj  =  /3j  and,  as  explained  in 
Section  2,  the  last  row  of  the  Toeplitz  matrix  is  the  reverse  of  the  first  one,  so  its  computation  may 
be  omitted,  hence  reducing  the  operation  count  by  fifty  percent.  The  Schur  algorithm  for  this  case 
follows: 

Schur  Algorithm  for  Symmetric  Positive-Definite  Matrices 


(C  *<0)  ...  &)=(*  it  ...  «n-l) 
1  <  fc  <  n  —  1,  Pfc  =  «(fc“*+1)/to*-I). 
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CA,  \ 


4.  The  Hyperbolic  Cholesky  Algorithm 

The  hyperbolic  Cholesky  algorithm  [4]  (see  also  [1,  7])  computes  both  Cholesky  factorisations 
A  =  UUT  and  A  =  CCT ,  where  UT  =  D\J*Ut  and  CT  =  D]j2LT,  of  the  symmetric  positive-definite 
matrix  A.  If  D  is  a  diagonal  matrix  whose  diagonal  equals  the  diagonal  of  A,  then 


4.1.  Connection  to  the  Generalised  Bareiss  Algorithm 


Premultiplying  both  sides  of  the  equation  by 


D 


1/2 


D 


1/2 


of  the  factorisation,  the  following  relations  between  the  quantities 
Hyperbolic  Cholesky  algorithms: 


gives,  due  to  the  uniqueness 
computed  by  the  Bareiss  and 


Thus,  as  also  shown  in  [5],  the  Hyperbolic  Cholesky  algorithm  is  a  ‘scaled’  version  of  the  generalised 
Bareiss  algorithm. 

Instead  of  the  ‘rotations’  of  the  generalised  Bareiss  algorithm,  the  Hyperbolic  Cholesky  algo¬ 
rithm  uses  hyperbolic  rotations  to  eliminate  a  pair  of  elements.  In  particular,  if  the  generalised 
Bareiss  algorithm  removes  element  (i  +  k,  i)  in  the  upper  matrix  and  element  (i,  i  +  k)  in  the  lower 
matrix  by  applying 

/  1  -c*i,i+k  \ 

\  ~&i+k,i  1  / 


9 


to  rows  i  +  k  and  i,  then  the  Hyperbolic  Cholesky  algorithm  removes  elements  in  the  same  positions 
via  the  hyperbolic  rotation 


1  (  1  ~Pi+k,i  \ 

y/l  ~  Phk,i  '  ~Pi+kli  1  ' 


where  pt+k,i  =  y/cti,i+kPi+k,i •  This  can  easily  be  seen  in  a  2  X  2  example,  where 


■(:  :)■ 


The  generalised  Bareiss  algorithm  computes 
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b  c 
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\b  c  ) 


(  a 


0  c  —  b2/a 

a  —  62/c  0 

b  c 


where  a  =  b/a  and  (3  =  bjc.  The  Hyperbolic  Cholesky  algorithm  computes 


f  1 
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where  p  =  b/y/ac.  Thus,  p  =  y/af]. 

The  process  of  forward  elimination  and  backsubstitution  proceeds  as  in  the  generalised  Bareiss 
but  with  Q  replaced  by 

Q"(D  £-1/2)' 

As  a  consequence,  the  Hyperbolic  Cholesky  algorithm  gives  rise  to  a  scaled  Schur  algorithm  based 
on  rotations 


1  (  1  ~P*\ 


4.2.  View  of  Hyperbolic  Cholesky  Factorisation  as  a  Downdating  Process 

Unlike  the  generalised  Bareiss  algorithm  the  Hyperbolic  Cholesky  algorithm  can  take  advantage 
of  the  symmetry  of  non-Toeplitz  matrices.  If  only  one  Cholesky  factor  is  desired,  it  suffices  to 
compute 

-r2 

where  A+  consists  of  the  upper  triangle  of  A  including  the  diagonal,  and  A#  consists  of  the  strict 
upper  triangle  of  A  excluding  the  diagonal.  As  a  product  of  hyperbolic  rotations,  the  matrix  Q h 
is  pseudo-orthogonal,  that  is. 


Ql 


(7  -/)«H=(7  _,) 
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(b)  Removing  Successive  Rows. 

Figure  4:  Two  Different  Orders  of  Rotations  in  the  Hyperbolic  Cholesky 
Algorithm  Applied  to  a  4  X  4  Matrix. 


Multiplying  the  previous  equation  by  its  transpose  and  exploiting  the  pseudo-orthogonality  of  Q h 
yields  CCT  =  VVT  -  WWT,  where  V  =  D~l'2A+  is  an  upper  triangular  matrix  with  positive 
diagonal  elements,  and  W  =  D~1/2A#  is  a  strictly  upper  triangular  matrix.  Denote  by  /,  the 
columns  of  W,  so  CCT  =  VVT  —  l, if ;  and  V(°)  =  V.  We  can  view  the  Hyperbolic  Cholesky 
algorithm  as  computing  a  sequence  of  Cholesky  downdating  problems:  each  problem  consists  of 

'  y(«'+i) 


premultiplying 


(7) 


by  hyperbolic  rotations  resulting  in 


(T) 


and 


y(«+i)y(«+i)r  _  y(')y(')T  _  [.[T 

such  that  y(‘+1)  is  an  upper  triangular  matrix  with  positive  diagonal  elements,  see  for  example  [11]. 
The  Hyperbolic  Cholesky  algorithm  computes  the  rotations  in  the  same  order  as  the  generalised 
Bareiss  algorithm.  This  is  also  the  order  that  a  downdating  process  would  choose,  as  Figure  4 
illustrates.  Instead  of  removing  successive  diagonals,  as  shown  in  Figure  4(a),  one  can  also  pipeline 
the  rotations  so  as  to  remove  successive  rows;  see  Figure  4(b). 

Reducing  the  Hyperbolic  Cholesky  algorithm  to  a  sequence  of  downdating  problems  may  facil¬ 
itate  its  round-off  error  analysis  due  to  the  availability  of  round-off  error  analyses  for  downdating 
problems. 


5.  The  Generalised  Levinson  Algorithm 

The  first  explicit  algorithm  for  solving  n  x  n  Toeplitz  systems  with  0(n2)  operations  was 
introduced  by  Levinson  in  1947  [12].  More  recently,  Levinson’s  algorithm  was  extended  to  non- 
Toeplitz  matrices  by  Delsarte,  Genin  and  Kamp  [6],  just  as  Bareiss  algorithm  [2]  was  extended  to  the 
generalised  Bareiss  algorithm  [5].  In  previous  work  [10]  we  have  shown  that  the  Schur  algorithm  can 
be  derived  as  the  result  of  trying  to  increase  the  degree  of  parallelism  in  the  Levinson  algorithm 
(i.e.  by  getting  rid  of  the  inner  products).  Conversely,  we  will  now  show  that  the  generalised 
Levinson  algorithm  is  a  simple  consequence  of  the  generalised  Bareiss  algorithm.  Our  version  of 
the  generalised  Levinson  algorithm  is  simpler  and  more  intuitive  than  the  one  in  [6]  as  it  does 
without  the  exchange  matrix  J,  and  it  establishes  a  direct  connection  to  the  other  factorisation 
methods. 
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5.1.  From  the  Generalised  Bareiss  to  the  Generalised  Levinson  Algorithm 
The  generalised  Bareiss  algorithm  computes 


Q 


CMS)- 


where  A  =  LU  and  A  =  UL  are  the  triangular  factorisations  of  a  non-singular  matrix  A.  Postmul- 
tiplying  both  sides  of  the  equation  by  ^  ^  ^  gives 

•CMfr'MK) 

since  A-1  =  U~l L~l  and  A~l  =  L~lU~l.  Thus,  by  applying  the  matrix  Q  to  the  identity  matrix 
the  generalised  Levinson  algorithm  determines  explicitly  the  unit  triangular  factors  L~x  and  U~x 
of  A-1  and  implicitly  (by  forming  U  =  £-1A  and  L  =  U~1A )  the  non-unit  triangular  factors  of  A. 

Now  we  consider  how  the  elements  of  Q  are  actually  determined  by  the  generalised  Levinson 
algorithm.  If  initially  aj0^  =  a,-,  1  <  *  <  n,  where  a,  is  the  t'th  row  of  A  then  the  fcth  step  of  the 
generalised  Bareiss  algorithm  combines  the  (i  +  fc)th  row  of  the  upper  matrix  with  the  ith 

row  of  the  lower  matrix: 


1  <  i  <  n  —  k, 


where 


-  Jk-i)  ,J-k+i) 


Pt+k,t  -  ai,i+k  /ai+k,i+k 


to  remove  element  (i+h,  i)  in  the  upper  matrix  and  element  (*,  i+k)  in  the  lower  matrix.  Eventually, 
a[k  ^  is  the  kth.  row  of  U  and  a*~n+^  is  the  kth  row  of  L. 

This  implies  that  the  fcth  row  of  U A-1  =  L~l  is  ip[k A-1  and  the  fcth  row  of 

LA'1  =  I/-1  is  ip{~n+k)  =  a[~n+k^  A-1.  In  particular,  =  aiA-1  =  ex  and  =  anA_1  =  ej, 
where  e,  are  the  n  x  1  canonical  vectors  with  a  one  in  position  t  and  zeros  everywhere  else.  This 
motivates 

VT  =  m_1  = «?', 

if 

0,-+*  =  ai+iA_1’  ip{~k)  =  a\~k)  A~l ,  1  <i<n-k, 

then  it  follows  by  induction  that 


',1°)  =  n  .  a  1  —  M  1  <  t  <  n. 


oti'i+k  =  Tl>i+kl)Aei/4’i  fc+1) Ae,-,  A+i fc,«  =  k+l)  Aei+klxp\k+kx)  Aei+k, 


and 


t{,+k  )  1  ~ Oi, i+k  ^  (  WJ 

tk)J  "v-ft+w  1  Aw- 


1  <  t  <  n  —  k. 


Note  that  only  elements  i  ...i-1-  k  of  0-*^  and  ip\  k '  may  be  non-zero.  To  keep  the  presentation 
simple,  we  do  not  take  advantage  of  this  fact  in  the  description  of  the  generalised  Levinson  algorithm 
below. 
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The  above  derivation  shows  that,  in  contrast  to  the  presentation  of  the  generalised  Levinson 
algorithm  in  [6],  there  is  no  need  for  the  exchange  matrix  J.  Its  occurrence  in  the  symmetric 
Levinson  algorithm  is  an  artefact  of  the  symmetric  Toeplitz  structure,  as  explained  in  the  next 
section. 

We  further  observe  that  the  inner  products  in  the  denominators  of  «,.,+*  and  /3,+i,,  can  be 
replaced  by  recursive  equations,  thus  saving  about  0(2k)  operations  in  step  k  of  the  generalised 
Levinson  algorithm.  Define  the  denominators  of  a, •,,•+*  and  /?;+&,{  by 

dtM'  =  ^W)Ae„  <£;■>  S  ^Aei+h, 


so  that 

<*M+fc  =  &'+*.«  =  i’i~k+1) • 

We  can  now  derive  recursive  equations  for  d\~k^  and  in  terms  of  d\~k+1^  and  initially 

d-0)  =  =  ej Ae;  =  a,-,-,  1  <  i  <  n. 

/  L\ 

From  the  above  definition  and  the  computation  of  xl>\  '  it  follows  that 

4-‘>  -  *!'*'-»«(  -  =  4-»"  -  <wW|-w) 

Similary,  one  shows  that 

d\+lc  =  (I  ~  ai,i+kPi+k,i)d\+k  \ 

Our  formulation  of  the  generalised  Levinson  algorithm  is  summarised  below. 


Factorisation  in  the  Generalised  Levinson  Algorithm 


1  <  t  <  n,  t/^0)  =  ej ,  d[0)  =  an 

1  <  k  <  n  -  1,  1  <  i  <  n  -  k, 


The  kth  rows  of  L~x  and  U~x  are  respectively  given  by  rp\.k  ^  and  n+k\  and  the  fcth  diagonal 
elements  of  U  and  L  by  d^k~1^  and  d^~n+k\ 

Because  it  computes  the  same  quantities  a,,,+*  and  /?,• hence  the  same  matrix  Q ,  as  the 
generalised  Bareiss  algorithm,  the  generalised  Levinson  algorithm  does  not  break  down  as  long  as 
all  contiguous  principal  submatrices  are  non-singular.  Moreover,  forward  elimination  and  backsub- 
stitution  can  be  done  in  the  same  manner  as  for  the  generalised  Bareiss  algorithm. 

For  a  symmetric  positive-definite  matrix  A ,  the  generalised  Levinson  algorithm  determines 
L_1,  U~x,  Di  and  D\j.  Alternatively,  the  Hyperbolic  Cholesky  algorithm  can  be  used  to  derive  a 
symmetric  version  of  the  generalised  Levinson  algorithm  that  computes  the  Cholesky  factorisations 
A-1  =  irTU~x  =  C~T  C~x  of  the  inverse. 
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5.2.  Specialisation  to  Toeplitz  Matrices 

If  the  matrix  T  at  hand  is  a  Toeplitz  matrix  then  the  generalised  Levinson  algorithm  reduces 
to  the  Levinson  algorithm  [12].  The  Levinson  algorithm  is  derived  by  applying  to  the  Bareiss 
algorithm  the  relations 


. /,(*)  _a(*)  A-l 

Vi+k  "  ai+kA  ' 


1  <  i  <  n  —  k. 


Hence,  only  the  vectors  and  are  needed;  see  Section  3.4.  Since  the 

generalised  Levinson  algorithm  computes  the  same  multipliers  a,-,i+jfc  and  /?,+*,«  as  the  generalised 
Bareiss  algorithm,  the  Toeplitz  case  simplifies  to 

=  ai<i+k  =  4>(kk~11)Tel/rp[--k+1)Teu  pk  =  Pi+k  =  4>[-k+1)Tek+l/tp[kk+-1)Tek+1 . 

In  contrast  to  the  Bareiss  algorithm,  which  requires  the  three  different  elements  ^  and 

only  one  set  of  quantities,  and  i’l  for  instance,  suffices  to  compute  both  mul¬ 

tipliers  for  the  Levinson  algorithm.  Hence,  about  half  of  the  arithmetic  operations  in  the  Levinson 
algorithm  can  be  saved  when  only  one  of  the  factorisations  is  desired  (that  is,  the  computation 
of  the  vectors  with  subscript  n  in  the  algorithm  below  can  be  omitted).  For  simplicity,  we  write 
<f>[k~1]  =  and  <t>(nk+l)  =  and  we  also  exploit  the  fact  that  only  elements  i,  . . .,  i  +  k 

of  and  tl>\~k^  may  be  non-zero. 

Factorisation  in  the  Levinson  Algorithm 


1  <  t  <  n,  =  <t>^  =  1?  t o 

!<*<»  —  !, 


0 

Jr 

II 

I—* - — ' 

Oi. 

i  . . 

■  t.k 

f/dk, 

pk  =  v>(~ 

dk+ 1  =  (1  -  (*kPk)dk 

(  4‘>  \  _  t 

i 

> 

-otk 

\  (  0 

h—  x 

1  ' 
Jr 

1 

-Pk 

1  , 

)  Ui-fc+1) 

o  ; 

II 

^3 

i 

-Pk 

-ak  \ 
1  ) 

1 

0  ) 

tk  )T  / dk 


The  A;th  rows  of  L~l  and  U~x  are  respectively  given  by 

0„_fc )  and  (Ofc_i  4>n  , 

and  the  fcth  diagonal  elements  of  U  and  L  by  dk  and  dn-k+\. 

If  the  matrix  T  is  symmetric  positive-definite  then  i\~k^  =  and  -  4>n  and 

because  of  L=JUJ  also  <pi~k^  =  <f>\ k^J.  Hence  the  algorithm  simplifies  as  follows. 

Factorisation  in  the  Levinson  Algorithm  for  Symmetric  Positive-Definite  Matrices 


1  <  i  <  n,  <^i°*  =  1,  d\  =  t0 
1  <  k  <  n  -  1, 

Pk  =  4>  ifc_1,(<i  •••  tk)T /dk,  dk+i  =  (l  -  pk)dk 

l  n 


^'  =  (i  -n)^ 


l)J  0 
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6.  A  Parallel  Algorithm  for  Persymmetric  Systems 


During  the  usual  process  of  solving  a  linear  system  Ax  =  b  as  described  in  Section  2.1,  the 
factorisation  and  forward  elimination  part  can  be  performed  simultaneously.  The  backsubstitution 
part,  however,  cannot  be  overlapped  with  the  other  two  because  the  first  step  in  the  solution  of  the 
upper  triangular  system  requires  the  last  element  obtained  from  the  solution  of  the  lower  triangular 
system  in  the  forward  elimination. 

We  will  now  show  how  one  can  combine  the  Bareiss  and  Levinson  algorithms  in  order  to 
perform  all  three  steps  of  solving  a  linear  system  Ax  =  b  in  parallel,  where  A  is  a  persymmetric 
matrix.  It  appears  that  this  algorithm  when  applied  to  Toeplitz  matrices  is  identical  to  one  given 
in  [8],  but  we  give  a  much  simpler  description  here. 

Suppose  that  there  are  enough  processors,  i.e.  0(n2),  available;  and  that  a  division,  or  a 
multiplication  followed  by  an  addition  constitutes  one  arithmetic  operation.  The  computation 


Q 


(T  I  b\  _  (U 

\T  I  V'U 


L~x 

u -x 


)•  (:)-(£*.)• 


where  the  multipliers  c*k  and  0k  are  computed  as  in  the  Bareiss  algorithm,  performs  the  factorisation 
and  forward  elimination.  It  requires  2(n  —  1)  parallel  operations  because  the  computation  of  the 
multipliers  cannot  be  overlapped  with  their  application. 

Because  T  =  UD\jL\j  ~  LD^Ui  and  because  of  persymmetry  we  actually  have 

(U  L~x  c\_(DhUL  L~x  c\_(  DlUl  L~x  c\ 

\L  U~x  d)~\DvLv  U-1  d)~\JDLLTJ  JU£tJ  d)' 

and  x  =  U~lc  =  U^x  D^x  c.  Because  the  Levinson  algorithm  delivers  the  elements  of  D^,  the 
divisions  y  =  D^xc  can  be  performed  concurrently  with  the  forward  elimination  (possibly  with  a 
lag  of  0(1)  steps).  Note  that  the  elements  of  the  vectors  c  and  y  are  computed  in  the  order  1, . . . ,  n, 
and  the  rows  of  JU^TJ  in  the  order  n,  ...,1.  Consequently,  the  rows  of  U^TJ  =  (JU^X)T  are 
available  in  the  order  l,...,n,  so  the  columns  of  JU^X  are  available  in  the  order  l,...,n  (the 
premultiplying  matrix  J  just  permutes  the  rows  and  has  no  effect  on  the  columns).  Thus,  the  ith 
column  of  U£x  is  available  by  the  time  the  »th  element  of  y  has  been  computed  so  that  the  linear 
combination  x  =  <7£" ly  of  the  columns  of  1  can  be  performed  concurrently  with  the  factorisation 
and  forward  elimination  in  2n  +  0(1)  parallel  operations. 

The  algorithm  for  Toeplitz  matrices  can  be  implemented  on  0(n)  processors  [8]  in  2n  parallel 
operations,  but  seems  to  require  a  considerable  amount  of  broadcasting.  If  a  truly  systolic  imple¬ 
mentation  without  broadcasting  is  desired  then  we  believe  that  this  algorithm  has  the  same  time 
complexity  as  the  one  in  [4],  which  uses  the  matrix  Q  to  perform  backsubstitution,  see  Section  3.3. 
At  last  note  that  the  asymptotic  time  complexity  is  the  same  for  persymmetric  and  for  Toeplitz 
systems. 
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